function [out] = mcsim(in, cfg) % TODO: 与Testing Framework融合，设计MCTestCase子类
%MCSIM 蒙特卡罗仿真
%
% Tips
% # 本函数将蒙特卡罗仿真分为多组，每组包含多次独立的仿真
% # 每组蒙特卡罗仿真采用相同的参数子集（从全体参数中选取）及相同的参数分布类型
% # 每次仿真首先运行仿真函数，再运行仿真后处理函数，各次仿真的输出数据可保存在独立的MAT文件中
% # 所有仿真完成后运行蒙特卡罗仿真后处理函数，可在该函数中读取MAT文件并做后处理
% # 调用本函数前，建议将工作目录设为仿真函数所在目录
% # 在进行仿真时，本函数将改变工作目录，若需访问调用本函数时的工作目录，参见cfg.fcnSim的帮助
% # 本函数自动将后处理函数输出保存在工作目录的mcsimout.mat中
%
% Input Arguments
% # in.grpParIdx: 长度为仿真组数的元胞向量，每个元素为向量，分别为每组蒙特卡罗仿真包含的参数子集在全体参数中的序号，未包含在该子集中的参数为0
% # in.grpParDistType: 长度为仿真组数的元胞向量，每个元素为长度与in.grpParIdx对应元素相等的向量，分别为对应参数的分布类型，1表示常量，-1表示负常量，2表示正态分布，3表示单边均匀分布，-3表示对称均匀分布，其它值默认为常量
% # in.grpSimNum: 长度为仿真组数的向量，分别为每组蒙特卡罗仿真的仿真次数
% # in.par: 表格，包含name、refVal、unit三个变量，说明如下
% 	in.par.name: 长度为仿真的全体参数个数的元胞向量，每个元素为字符串，分别为各参数的名称
%   in.par.refVal:
%   行数为仿真的全体参数个数，列数为3的矩阵，各参数的参考值，第1/2/3列分别为当该参数分布类型为（负）常量/正态分布/均匀分布时，该参数的（负）全值/标准差/最大值，对应的分布律为（-）RefVal/N(0, RefVal)/单边(0, RefVal)或对称(-RefVal, RefVal)
% 	in.par.unit: 长度为仿真的全体参数个数的元胞向量，每个元素为字符串，分别为各参数的单位
% # cfg.fcnSim: 单次仿真函数句柄，该函数接收一个结构体输入参数，包含curGrpIdx、curGrpSimNum、curSimIdx、curSimParVal、curSimParDesc、parRefVal、CWD等域，其中：
%   curGrpIdx为标量，是本组仿真的组序号（从1开始）；curGrpSimNum为标量，是本组仿真包含的仿真次数；curSimIdx为标量，是本次仿真在本组内的次序号（从1开始）；
%   curSimParVal为长度为全体参数个数的列向量，是本次仿真的参数值，由本函数根据各参数的参考值和分布类型随机生成；curSimParDesc为字符串，包含本次仿真参数子集中各参数的信息；
%   parRefVal同in.par.refVal；CWD为字符串，是调用本函数时的工作目录
% # cfg.fcnPostProc: 蒙特卡罗仿真后处理函数句柄，在所有仿真完成后执行，该函数接收一个结构体输入参数，包含in的所有域，并有一个输出参数
%
% Input Arguments
% # out: cfg.fcnPostProc函数的输出参数
tmp_parTypeDesc = {'（minus）whole value', 'std', 'max'};
% 按组运行仿真
disp('*****************************************Started Monte Carlo Simulation**************************************************');
tic;
tmp_CWD = pwd;
addpath(tmp_CWD)
mid_cumSimNum = cumsum(in.grpSimNum);
parfor mid_curTotalSimIdx = 1:mid_cumSimNum(end) % 使各组各次仿真均能并行运行
    mid_curGrpIdx = findendingcompare(mid_cumSimNum, mid_curTotalSimIdx, @ge);
    tmp = [0; mid_cumSimNum(:)];
    mid_curSimIdx = mid_curTotalSimIdx - tmp(mid_curGrpIdx);
    % 生成本组仿真的参数说明
    mid_curGrpSimNum = in.grpSimNum(mid_curGrpIdx); %#ok<*PFBNS>
    mid_curSimParDesc = '';
    tmp_curGrpParIdx = in.grpParIdx{mid_curGrpIdx};
    tmp_curGrpParDistType = in.grpParDistType{mid_curGrpIdx};
    if length(tmp_curGrpParIdx) ~= length(tmp_curGrpParDistType)
        warning('mcsim:lengthmismatch', ['The length of grpParIdx and grpParDistType of group ' int2str(mid_curGrpIdx) ' are not the same.']);
    end
    for tmp_i=1:length(tmp_curGrpParIdx)
        tmp_parIdx = tmp_curGrpParIdx(tmp_i);
        tmp_parDistType = tmp_curGrpParDistType(tmp_i);
        mid_curSimParDesc = strcat(mid_curSimParDesc, [in.par.name{tmp_parIdx} '：' ...
            num2str(in.par.refVal(tmp_parIdx, abs(tmp_parDistType))) ...
            '（' in.par.unit{tmp_parIdx} '，' tmp_parTypeDesc{1, abs(tmp_parDistType)} '） ']);
    end
    % 按次运行仿真
    if isunix
        tmp_taskPath = [tmp_CWD '/task' int2str(get(getCurrentTask(), 'ID'))];
    else
        tmp_taskPath = [tmp_CWD '\task' int2str(get(getCurrentTask(), 'ID'))];
    end
    if ~isfolder(tmp_taskPath)
        mkdir(tmp_taskPath);
    end
    cd(tmp_taskPath)
    disp(['**************************************************Started the ' num2str(mid_curSimIdx) '-th Monte Carlo Simulation of Group ' num2str(mid_curGrpIdx) '**************************************************']);
    % 生成本次仿真的参数
    mid_curSimParVal = zeros(size(in.par.refVal, 1), 1);
    for tmp_i=1:length(tmp_curGrpParIdx)
        tmp_parIdx = tmp_curGrpParIdx(tmp_i);
        tmp_parDistType = tmp_curGrpParDistType(tmp_i);
        tmp_parRefVal = in.par.refVal(tmp_parIdx, abs(tmp_parDistType));
        switch(tmp_parDistType)
            case 2 % 正态分布
                mid_curSimParVal(tmp_parIdx, 1) = tmp_parRefVal * randn;
            case 3 % 单边均匀分布
                mid_curSimParVal(tmp_parIdx, 1) = tmp_parRefVal * rand;
            case -3 % 对称均匀分布
                mid_curSimParVal(tmp_parIdx, 1) = tmp_parRefVal * (rand-0.5) * 2;
            case -1 % 负全值
                mid_curSimParVal(tmp_parIdx, 1) = -tmp_parRefVal;
            otherwise % 全值
                mid_curSimParVal(tmp_parIdx, 1) = tmp_parRefVal;
        end
    end
    % 运行本次仿真
    cfg.fcnSim(struct('curGrpIdx', mid_curGrpIdx, 'curGrpSimNum', mid_curGrpSimNum, ...
        'curSimIdx', mid_curSimIdx, 'curSimParVal', mid_curSimParVal, 'curSimParDesc', mid_curSimParDesc, ...
        'parRefVal', in.par.refVal, 'CWD', tmp_CWD));
    disp(['**************************************************Finished the ' num2str(mid_curSimIdx) '-th Monte Carlo Simulation of Group ' num2str(mid_curGrpIdx) '**************************************************']);
    cd(tmp_CWD)
end
rmpath(tmp_CWD)
disp('**************************************************Finished All Groups of Monte Carlo Simulation**************************************************');
toc;
% 运行蒙特卡罗仿真后处理
close all
out = cfg.fcnPostProc(in);
save('mcsimout', 'out');
disp('**************************************************Finished Monte Carlo Simulation PostProc**************************************************');
toc;
end